clc
clear all


n = 150; % number of target group elements
ob_p = 0.2;   % observation probability
q = 1; % probability of not being corrupted by (multiplicative) outlier noise, q=1 means no outlier noise
delta = 1;    % noise level in additive noise, delta=0 means no additive noise

Number_of_d = 5;   % number of points
Array_d = 0 + 5 * (1:Number_of_d)';   % number of permutations 


Number_of_Sim = 3;   % number of repeats

% error 
err_SpeIni = zeros(Number_of_Sim, Number_of_d);
err_QR = zeros(Number_of_Sim, Number_of_d);
err_GPM = zeros(Number_of_Sim, Number_of_d);

% running time
time_Spectral = zeros(Number_of_Sim, Number_of_d);
time_QR = zeros(Number_of_Sim, Number_of_d);
time_GPM = zeros(Number_of_Sim, Number_of_d);

for index_d = 1:Number_of_d
    d = Array_d(index_d);    
    
    for index_sim = 1:Number_of_Sim
        fprintf('index_n: %d of %d, index_sim: %d of %d\n', index_d, Number_of_d, index_sim, Number_of_Sim);
        [C, row_num, P_true] = GenerateData_Perm(n, d, delta, ob_p, q);
        
        [t_Spectral, P_Spectral] = Spectral_Perm( C, d);
        [t_QR, P_QR] = QR_Perm( C, P_Spectral, row_num );
        [t_GPM, P_GPM] = GPM_Perm( C, P_Spectral );
        
        err_SpeIni(index_sim, index_d) = Error_Perm(P_true, P_Spectral);
        err_QR(index_sim, index_d) = Error_Perm(P_true, P_QR);
        err_GPM(index_sim, index_d) = Error_Perm(P_true, P_GPM);
                
        time_Spectral(index_sim, index_d) = t_Spectral;
        time_QR(index_sim, index_d) =  t_QR + t_Spectral;
        time_GPM(index_sim, index_d) = t_GPM + t_Spectral;
    end
end

% 
figure;
set(gcf,'units','points','position',[10,10,800,400])
% title('d=5, observe p_1 = 0.1, noise p_2=0.7')

subplot(1,2,1);
set(gca,'fontsize',15)
plot(Array_d,1-sum(err_SpeIni)/Number_of_Sim,'-kp','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_d,1-sum(err_QR)/Number_of_Sim,'-g*','LineWidth',2,'MarkerSize',10);
plot(Array_d,1-sum(err_GPM)/Number_of_Sim,'-ro','LineWidth',2,'MarkerSize',10);
hold off
legend('Eigen','QR','GPM','Location','northwest')
xlabel('$d$','Interpreter','latex')
ylabel('\textbf{Recovery Rate}','Interpreter','latex')


subplot(1,2,2);
set(gca,'fontsize',15)
plot(Array_d,sum(time_Spectral)/(Number_of_Sim),'-kp','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_d,sum(time_QR)/(Number_of_Sim),'-g*','LineWidth',2,'MarkerSize',10);
plot(Array_d,sum(time_GPM)/(Number_of_Sim),'-ro','LineWidth',2,'MarkerSize',10);
hold off
legend('Eigen','QR','GPM','Location','northwest')
xlabel('$d$','Interpreter','latex')
ylabel('\textbf{Time (s)}','Interpreter','latex')

sgtitle(['$n$ = ', num2str(n), ', $p$ = ', num2str(ob_p), ', $q$ = ', num2str(q), ', $\delta$ = ', num2str(delta)],'Interpreter','latex');